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We show how a fixed point based boundary-layer analysis technique can be used to obtain the 
steady-state particle density profiles of driven exclusion processes on two-lane systems with open 
boundaries. We have considered two distinct two-lane systems. In the first, particles hop on the 
lanes in one direction obeying exclusion principle and there is no exchange of particles between the 
lanes. The hopping on one lane is affected by the particle occupancies on the other, which thereby 
introduces an indirect interaction among the lanes. Through a phase plane analysis of the boundary 
layer equation, we show why the bulk density undergoes a sharp change as the interaction between 
the lanes is increased. The second system involves one lane with driven exclusion process and the 
other with biased diffusion of particles. In contrast to the previous model, here there is a direct 
interaction between the lanes due to particle exchange between them. In this model, we have looked 
at two possible scenarios with constant (flat) and non-constant bulk profiles. The fixed point based 
boundary layer method provides a new perspective on several aspects including those related to 
maximal/minimal current phases, possibilities of shocks under very restricted boundary conditions 
for the flat profile but over a wide range of boundary conditions for the non-constant profile. 



I. INTRODUCTION 

Many unique non-equilibrium phenomena such as 
boundary- induced phase transitions [l|, spontaneous 
symmetry breaking 0], phase separation [3| are ob- 
served in one-dimensional driven exclusion processes. 
These features that are very special to systems far-from- 
cquilibrium led to extensive studies of variety of driven 
processes. Major developments in this direction started 
with the asymmetric simple exclusion process (ASEP) in 
which particles, after being injected at one boundary at 
a specific rate, hop in a specific direction on a finite one- 
dimensional lattice obeying mutual exclusion rule d, Q . 
A steady particle current is sustained through withdrawal 
of particles from the other boundary at a given rate. This 
simple model was later followed by more general systems 
involving more than one species of particles > more 

complicated lattice [§l-fl4j or dynamics of particles [l5[ 
etc. All these models have a few common, basic features 
such as hopping of particles with bias in one direction and 
hence a nonzero particle flux, mutual exclusion between 
the particles and open boundaries with particle injection 
and withdrawal at given rates. In contrast to equilib- 
rium systems, these one-dimensional, driven, many par- 
ticle systems exhibit boundary induced phase transitions 
in their steady-state. In various phases, density profiles 
have distinct shapes which, for a given process, are com- 
pletely dependent on the boundary rates. The boundary 
rates, therefore, are the most natural variables for the 
phase-diagram representing the phases and phase transi- 
tions. 

In order to characterize the phase transitions, it often 
appears convenient to look at the steady-state density 
profile which describes the average particle occupancy 
of various sites. Density profiles have interesting shapes 
with extended bulk parts and one or more boundary layer 
parts which are narrow regions over which the density 
varies rapidly. The nature and location of the bulk, as 



well as the boundary layer parts of the profile, change as 
the boundary rates change. For example, the boundary 
layers may be located near one or both the boundaries or 
may appear in the interior of the lattice. The latter one, 
commonly known as a shock, separates high and low den- 
sity bulk regions [l6j . That the boundary layer need not 
be confined to the boundary alone and the formation of a 
shock can be well characterized through a deconfmement 
transition of the boundary layer from the boundary have 
been shown in [l7l - IT9| . In the next section, we shall pro- 
vide a more technical description of the boundary layer 
in terms of an appropriate differential equation. 

Recent studies show that the information about 

the bulk density can be obtained by studying the bound- 
ary layer parts of the density profile. Although methods 
of boundary layer analysis [20( allow us to find the ana- 
lytical expressions of the boundary layers [l7[, in many 
cases, it becomes technically challenging to obtain an- 
alytical expressions for the boundary layers. However, 
since the boundary layers are expected to merge to the 
bulk in the appropriate limit, one may obtain informa- 
tion about the bulk density by studying the fixed points 
of the differential equations describing the boundary lay- 
ers. For constant bulk profiles, this method is especially 
powerful since, in this case, the bulk densities must be 
the same as the fixed point values of the boundary layer 
equation. 

In the present work, we apply fixed point based bound- 
ary layer method Q to two driven many-particle systems, 
each composed of two one-dimensional lattices or lanes. 
Every lane has N sites and particles can hop to the neigh- 
boring sites on the lane obeying specific rules that are 
mentioned in detail in the respective sections. The par- 
ticle dynamics of the two models considered here differ 
due to the following reasons. In the first model, parti- 
cles do not hop from one lane to the other but the hop- 
ping rate on one lane changes based on the particle occu- 
pancy of the neighboring sites on the other lane 0, HI, [lj| . 
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In contrast, the second model involves two lanes which 
can mutually exchange particles. While particles on one 
lane undergo asymmetric simple exclusion process with 
hopping in a specific direction, those on the other lane 
undergo biased diffusion [l(| EH, EH . Such driven sys- 
tems have similarities with intracellular transport pro- 
cesses in which molecular motors, receiving energy from 
the hydrolysis of adenosine triphosphate (ATP), move 
in a particular direction along microtubules. The first 
model may mimic the situation where a molecular motor 
moving along one channel with a large cargo attached to 
it creates an obstruction for the motion of other molec- 
ular motors on the neighboring channel. A second lane 
with particle diffusion on it has been introduced earlier to 
represent the environment in which molecular motors dif- 
fuse during t he p eriod when they are not attached to the 
microtubule [10( ■ More recently, similar two-lane models 
have been introduced to describe extraction of membrane 
tubes by molecular motors [IH, EU • ^ n both the models 
that we study here, lanes are coupled to boundary reser- 
voirs that maintain specific densities at the two ends of 
a lane. Our primary aim is to illustrate how the method 
can be used to predict steady-state density profiles under 
different boundary conditions and to quantify various fea- 
tures of the boundary layers including the height of the 
shock, approach to the bulk etc. 

The motivation behind this study stems from the fact 
that till now there exists no general framework to study 
phase transitions of this small class of non-equilibrium 
systems. Although this is a much bigger issue concern- 
ing the entire subject of non-equilibrium statistical me- 
chanics, lack of a general framework even for these driven 
systems is surprising. Due to the presence of more than 
one density variable, driven multi-lane systems 
in general, are technically challenging. Previous stud- 
ies on two lane systems, with particle occupancy on one 
lane affecting the hopping on the other, show that this 
issue requires a generalization of the extremal current 
principle [l4[ which has been used earlier to predict the 
phase diagram of a single lane ASEP involving a sin- 
gle density variable [H, [2l[ . Inadequacy of mean-field 
analysis for certain cases has motivated development of 
cluster approximation Q and all these studies have re- 
vealed existence of new interesting phases including sym- 
metry breaking. Studying the stability properties of the 
bulk plateaux, a variety of phase diagrams has also been 
found for systems with ASEP on one lane coupled to dif- 
ferent kinds of particle's motion on the other lane fl3j . 
All these observations motivate us to verify the versatil- 
ity of the boundary layer based method to this class of 
problems. The present paper is restricted to only two 
specific models from this class. It would be interesting 
to apply this method to study phase diagrams of other 
variants of two lane systems. This will be discussed in 
our future publications. 

The present method appears to be general since it does 
not rely on any explicit analytical solution of the density 
profile, yet it provides a lot of physical insight regard- 



ing the location of the boundary layer, value of the bulk 
density, nature of the phase transitions etc; all these ob- 
tained analytically through a phase-plane analysis of the 
differential equation describing the boundary layer. We 
predict the shape of the entire density profile from the 
fixed points and the phase portraits of the boundary layer 
equation. Most importantly, we find that various known 
results, obtained through development of different meth- 
ods and hypotheses including extremal current princi- 
ple, follow as natural consequence of the phase-portrait 
of the boundary layer differential equation. There are 
resemblances between the present method and the ex- 
tremal current principle but it would require more work 
to establish a direct connection which might provide a 
natural basis for the hypothesis made for the extremal 
current principle. 

The stochastic dynamics of particles undergoing ASEP 
can be described through discrete master equations. For 
large ./V (N —> oo), and small lattice spacing, a (a — > 0) 
with Na finite, we may go over to a continuum limit in 
which a lattice site i is replaced by a continuous posi- 
tion variable x = i/N. In the continuum, long time, long 
length scale limit (the so called hydrodynamic limit), a 
statistically averaged master equation appears like a con- 
tinuity equation describing the time-evolution of a den- 
sity variable, say, p{x, t), in terms of the particle current 
j. Since bulk phase transitions are large length scale 
phenomena and only certain gross features are crucial 
at these length scales, a continuum formulation is ex- 
pected to be sufficient for our purpose. Our boundary- 
layer method is applicable to the steady-state version of 
such continuum equation. It is important to note that 
this continuum approach, however, produces a narrow 
boundary layer whose width varies with N. This will 
be shown explicitly in the next section by obtaining the 
boundary-layer solution for the simplest ASEP model in- 
troduced at the beginning. 

An explicit derivation of the steady-state hydrody- 
namic equation requires obtaining the stationary flux, 
j(p) from the microscopic dynamics of the model. For the 
first model, this derivation is simple due to spatially un- 
corrected nature of the steady-state and for the sec- 
ond model, we use a mean-field current-density relation. 
Since our primary aim is to elucidate how the fixed point 
based boundary-layer method works, we have used mean- 
filed approximation for convenience. It will be clear from 
the analysis that the method is robust in the sense that 
it works equally well for other hydrodynamic equations 
irrespective of the approximation that has been used to 
derive these. Hence, the outcome of this method is exact 
to the extent to which the starting hydrodynamic equa- 
tion is exact. A phase-plane analysis for the first model 
shows that, for certain boundary conditions, the density 
profile is strongly influenced by a saddle fixed point of 
the boundary layer differential equation and finally, as 
a consequence of this, the bulk profile changes drasti- 
cally to a new value as the interaction between the lanes 
is increased. The second model provides an interesting 
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fixed point diagram with saddlc-nodc bifurcations of the 
fixed points [22j of the boundary layer equation. These 
bifurcations appear at special densities that correspond 
to the maximum or minimum of the particle current and 
maximal/ minimal current phases appear as natural con- 
sequences of the flow trajectories towards these special 
densities. 

In section II, we discuss a few general properties of the 
boundary layer solution. Section III presents the phase 
plane analysis of the first model. In section IV, we con- 
sider the second model with two distinct cases of constant 
and non-constant bulk profiles. We summarize our main 
findings in section V. 

II. BOUNDARY LAYERS 

In order to illustrate some of the basic features of the 
boundary layer, here we choose the simplest model of 
ASEP in which particles hop in a particular direction on 
one lane obeying the exclusion rule. In addition, we also 
include processes that involve adsorption and evapora- 
tion of particles to and from the lane at rates proportional 
to uj a and uid, respectively. The lane is coupled to bound- 
ary reservoirs which maintain fixed particle densities, pi 
and p r , at left and right boundaries of the lane. For the 
following equation, equal particle adsorption and evapo- 
ration rates are assumed. In a continuum mean-field de- 
scription with the lattice size scaled to unity {Na = 1), 
the density variable t) satisfies the differential equa- 
tion 

|-.g + (*-i)g + n(i-„>. (D 

Here O = u> a N = uJdN and e = ™ is a small parameter. 
In the absence of adsorption and evaporation processes, 
the hydrodynamic equation with only hopping of parti- 
cles can be expressed in the form of a continuity equa- 
tion If = "^(~ e lf + 3 P ), wh erc j p = p(l - p) is the 
particle current associated with the hopping process. In 
addition, there is a diffusive current which appears along 
with a prefactor, e, which becomes small as N — > oo. In 
order to obtain the steady-state density profile, we need 
to solve 

4^ + (2p-l)^ + fXl-2p) = (2) 
ax * ax 

with boundary conditions p(x = 0) = p% and p(x — 1) = 

Pr- 
ill the limiting case, e = 0, equation @ becomes a 
first order equation and its solution cannot, in general, 
satisfy two boundary conditions. The vanishing higher- 
order derivative term, also known as the regularization 
term, helps avoiding the singularity in the differential 
equation. In these problems, boundary layers are ex- 
pected to appear near either one of the boundaries or in 
the interior of the lane. Based on the values of the bound- 
ary densities, different solutions arise and the techniques 



of the boundary-layer analysis allow one to obtain uni- 
form approximation to the solution of @ order by order 
in e under given boundary conditions. For very small e, 
the second order derivative term of equation ([5]) can be 
neglected. The solution of the resulting first order equa- 
tion describes the major part of the density profile. This 
solution, known as the outer solution in the boundary- 
layer language, is referred in the following as the bulk 
solution. For this particular example, the bulk solution 
is 

Po(x) = flx + C (3) 

with Co being the integration constant whose value can 
be determined from the boundary condition that this so- 
lution satisfies. For example, for a density profile with 
the bulk part satisfying the boundary condition at x = 0, 
Co = pi ■ It is possible that along with this bulk solution, 
a boundary layer appears near x = 1 (see figure [T] for 
a typical density profile that appears in the low density 
phase for the ASEP under consideration [l|| U3] ) • The 
boundary layer here satisfies the boundary condition at 
x = 1 and merges to the bulk solution at the other end. 
In order to satisfy two conditions, the second derivative 
term of ([2]) becomes necessary for the description of the 
boundary layer. Thus, in general, higher order deriva- 
tive terms dominate the behavior in the boundary layer. 
To focus on the boundary layer, we introduce a rescaled 
variable x = x ~^° ; where Xq, which is arbitrary at this 
stage, specifies the location of the boundary layer after 
appropriate boundary conditions are implemented. In 
terms of this rescaled variable, ([2]) appears as 

^ + ( 2/ >-l)^+eO(l-2p)=0. (4) 

aar ax 

The O(e ) solution of the boundary layer can be obtained 
from 

Since, the boundary layer is narrow, this also implies that 
at this order, the particle non-conserving processes have 
negligible influence on the boundary layer. Irrespective 
of their appearance in the interior or the boundary of 
the lane, solution of such equations are, in general, re- 
ferred as the boundary layer solutions in the following. 
If the boundary layer appears near x = 1, we expect the 
boundary layer solution to merge to the bulk at x — > — oo 
and satisfy the boundary condition p(x = 0) = p r . For 
boundary layers appearing in the interior of a lane sep- 
arating high and low density bulk solutions, we expect 
the boundary layer solution to merge to appropriate bulk 
densities in the x — > ±oo limits. 
Integrating ([5]) once, we have 

f_ +{p *-p) = co, (6) 

where cq is the integration constant. The saturation of 
the boundary layer to the bulk density, pi, = p Q (x = 1), 
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requires cq = p\ — p\>. This relation between cq and pb 
fixes the range of physically acceptable values of Cq. Ap- 
proach to different bulk values can be also understood by 
obtaining phase-portraits of such equations for different 
values of Co in this range. Fixed points, p* , of this equa- 
tion are the solutions of the equation p* 2 — p* — cq = 0. 
It often appears convenient to plot these fixed points as 
functions of cq. Obtaining the stability properties of the 
fixed points, one may indicate on this diagram, how, 
starting from a given initial density, the solution flows 
towards or away from a given fixed point. Such diagrams 
will be referred in the following as fixed point diagrams. 

Instead of elaborating on the fixed point diagram of 
this model in detail, we discuss some of the properties of 
the boundary layer by focusing on the density profile of 
figure ([T]). Clearly, the boundary layer appearing at x = 1 

0.5i ' ' ' ' 1 




x 



FIG. 1: A typical density profile in the simplest ASEP model 
considered in this section. 



goes beyond the physical size of the lane. In case the 
boundary layer has a stable fixed point, which is the case 
for the tanh type boundary layer shown in figure (fTj). 
one may expect to see a shock as the boundary condi- 
tion at x = 1 is appropriately adjusted. The boundary 
layer just deconfines from the x = 1 boundary and en- 
ters into the bulk in the form of a shock when the stable 
fixed point of the boundary layer is exactly the same as 
the boundary condition at x = 1. At this instance, the 
x — >• oo end of the boundary layer is just in the phys- 
ical region. Now, if the boundary condition at x = 1 
is raised to a slightly higher value, the boundary layer, 
that has already reached its saturation at x = 1 cannot 
anymore satisfy the boundary condition. It is, in this 
case, that the boundary-layer deconfines and enters into 
the bulk as a shock with the i — > oo merging to an- 
other high-density bulk solution. This example of ASEP 
presents a simpler problem since, in this case, one can ex- 
plicitly solve the boundary layer equations. Fixed point 
diagrams are especially useful for those cases for which 
such exact solutions cannot be obtained. 



III. MODEL-1: TWO-LANE SYSTEM 
WITHOUT PARTICLE EXCHANGE 

In this model, we have two lanes on which particles 
hop unidirectionally at rates as shown in figure 




has to merge to the bulk density, pb- We expect the bulk- 
density pb to be an unstable fixed point of the boundary 
layer to which the boundary layer merges as x — > — oo. 
Whether the other end of the boundary layer (x — > oo 
limit) approaches a fixed point or flows indefinitely is 
decided by the nature of the boundary-layer equation and 
can be seen clearly from the fixed point diagrams. In this 
particular case, boundary layer equation ([5]) has simple 
solutions 0, El 

Pu{x)= 1 -+ [1 ~^ Pb) tanh(^+0 and (7) 

Pu(x) = i+(I^coth(A+£), (8) 

where £ is a constant that gives the center of the bound- 
ary layer and w = 1 J > 2pb is the width of the boundary 
layer. While w describes the approach of the boundary 
layer to the bulk asymptote, £ helps visualizing the shock 
formation. These solutions also show how the boundary 
layers scale with the system size, N. The boundary layer 
of figure ([TJ is described by the tanh solution which satis- 
fies the boundary condition at x = 1 before approaching 
its stable fixed point, 1 — pi,, as x — ¥ oo, a limit that 



H h • h H • + * 

FIG. 2: The hopping rates on one lane depend on the total 
particle occupancy of the neighboring sites on the other lane. 
On a lane, a particle hops to the empty, forward site at rates 
1, (1 + /3)/2 and /3 if the neighboring sites on the other lane 
are empty, half-filled and completely filled respectively. 



Since the hydrodynamic equations have been already 
derived from the microscopic dynamics of the model, we 
quote these equations from the previous literature Q . In 
the steady-state, these equations are 

_^ + e i.[(l_ m(T) ^ ]= 0, and (9) 
ox ax ax 

where e = j-^ and p and a are the average densi- 
ties on the lanes with the corresponding currents j p = 
p(l - p)[l - (I - 0)a] and j a = ail - a)[l - (1 - (3)p\. 
Here, /3 measures the strength of the interaction between 
the two lanes with (3 = 1 representing the non-interacting 
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case. From now onwards, we use m = 1 — /3. These equa- 
tions are to be supplemented with the boundary condi- 
tions (p(x = 0) = pi, p(x = 1) = p r ) and (a(x = 0) = 07, 
a(x = 1) = a r ) at the two ends of the lanes. Instead of 
considering nonlinear regularization terms, we follow a 
phenomenological approach and choose simpler regular- 

ization terms of the form e^jf and e^gr- After studying 
the steady-state for this simpler situation, we argue that 
no new feature emerges with the actual regularization 
terms. 

Hydrodynamic equations that we study in the follow- 
ing are, 



Oj 



d 2 p 
OX ox z 

dx dx 2 



(11) 
(12) 



Equations ([TTj) and (|12|) admit constant solutions which 
correspond to constant bulk profiles. For the boundary- 
layer solutions, we re-express these equations in terms of 
x and integrate once to obtain 



dp 

— = p(l - p)(l -ma) +d, 

3(7 

— = a(l - cr)(l - mp) + d x . 



(13) 
(14) 



Here c\ and d\ are the two integration constants. The 
saturation of the boundary layers to the bulk densities 
Pb and (if, is ensured through the choice 



ci = -pb(l - Pb)(l - mcTb), and 
d\ = -<Tb(l - c&)(l - mpb). 



(15) 
(16) 



To begin with, we consider a situation where a density 
profile has a constant bulk part satisfying one boundary 
condition and a boundary-layer part satisfying the other 
boundary condition. As discussed earlier, one may con- 
clude that the bulk density value which is also a bound- 
ary density at one end, is a fixed point of the bound- 
ary layer equation. The boundary layer, in this case, is 
the solution of (fT3|) and (|T4| that starts from an initial 
density, which, here, is the boundary condition satisfied 
by the boundary layer, and approaches the fixed point. 



Hence, for a set of boundary conditions, we may find 
out all the fixed points of the boundary layer equations 
with ci and d\ fixed by assuming bulk densities to be 
same as one set of the boundary values. Whether this 
choice of bulk density and the corresponding boundary 
layer form an acceptable solution for the density profile, 
depends on the stability properties of fixed points and 
the phase-plane trajectories. Although settling this issue 
is much simple when the boundary layer equation has 
only stable or unstable fixed points, this is not so when 
the boundary layer equation has a saddle fixed point in 
addition to stable or unstable fixed points. In order to 
see these features, the most general approach would be 
to obtain the fixed point diagram as functions ci and d\ . 
For finding out the fixed points, p* and cr*, one has to 
solve algebraic equations, p*(l — p*)(l — ma*) + c\ = 
and a*(l — c*)(l — mp*) + d\ =0. Solving a fifth order 
polynomial equation for either p* or cr*, one may find 
five fixed points of which only real, positive fixed points 
of values p*(and a*) < 1 are the physically acceptable 
ones. Clearly, this would require a multi-dimensional pa- 
rameter space on which the topology of the manifolds 
describing the fixed-point would be displayed. Instead of 
this detailed fixed point diagram, one may also find fixed 
points of (fH?)) and (|T4"]) with values of Ci and c?i obtained 
for all possible combinations of pb and ab- This allows 
us to have a generic picture for the number of physically 
acceptable fixed-points and their stability properties on 
the entire pb — ab plane. 



Instead of providing a detailed picture, we consider a 
specific case of bulk densities pb — 0.4, and ab = 0.01. 
It turns out that this corresponds to the same shape of 
the density profile we have just discussed. Out of five 
sets of fixed points of the boundary layer equation, only 
two sets are physically meaningful and these fixed points 
govern the phase-plane trajectories of the boundary layer 
solutions. The other fixed points are not crucial for our 
analysis since these are either imaginary or are unphys- 
ical (fixed point values larger than 1). Table 1, in the 
following, provides a list of values of ci, d±, two physi- 
cally acceptable sets of fixed points (p\,a*) and (p*.; ^); 
and the corresponding sets of eigenvalues (An, A12) and 
(A21, A22) for different values of m. 



m 


Cl 


di 




(/>2,02) 


(An, A12) 


(A21, A22) 




0.48 
0.51 
0.58 
0.68 


-0.23884 
-0.23877 
-0.238608 
-0.238368 


-0.00799 
-0.00788 
-0.00760 
-0.0072072 


(0.4,0.01) 
(0.4,0.01) 
(0.4,0.01) 
(0.4,0.01) 


(0.5992,0.01134) 
(0.5991,0.0114) 
(0.598, 0.01178) 
(0.5980,0.0122) 


(0.7927,0.1981) 
(0.781,0.198) 
(0.754,0.197) 
(0.71556, 0.1965) 


((0.698,-0.199) 
(0.6808,-0.1997) 
(0.638,-0.195) 
(0.5794,-0.20) 





r 



The eigenvalues show that of the two physically mean- phase plane trajectories toward the unstable fixed point 
ingful fixed points, (pi, a*) is an unstable fixed point and (/?* = 0.4, cr* = 0.01) in the x —5- —00 limit for various 
(pgjCrj) is a saddle one. Figure [3] shows the flow of the initial values. Although, one can attempt to draw the 
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FIG. 3: Phase-plane diagram on the p — a plane for dif- 
ferent initial conditions with m = 0.4. Initial conditions 
(p(x = 0),a(x = 0)) as we move towards the outermost line 
are (0.2,0.4), (0.2,0.6), (0.2,0.7), (0.2,0.8). The arrows on 
the lines, indicate the direction of more negative values of x. 
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0.2 0.4 x 0.6 0.8 



FIG. 4: Density profiles (p) plotted for m — 0.4 and left 
boundary conditions (pi,ai) = (0.4,0.01). Right boundary 
conditions, here, are same as the initial conditions of figure 
131 The lower most profile corresponds to (p r ,a r ) — (0.2,0.4). 
The rest follow the same order as specified in figure[3j All the 
plots are obtained with e = 0.006. 




0.4 x 0.6 



FIG. 5: Density profiles, a, plotted for m = 0.4 and same 
boundary conditions as that of figure [4] The lower most pro- 
file corresponds to (p r ,a r ) — (0.2,0.4). The rest follow the 
same order as specified in figure [3] All the plots are obtained 
with e = 0.0062. 



(p,cr) — (0.2,0.8) comes to close proximity of the sad- 
dle fixed point. This is reflected in the boundary layer 
of the density profile plotted in figure [4] with boundary 
condition, (p r ,a r ) = (0.2,0.8). We next study how a 
trajectory with a given initial condition changes with 
m. Figure [5] shows phase plane trajectories with ini- 
tial condition (0.2,0.8) for different values of m. From 




flow trajectories intuitively, flow trajectories in the fig- 
ures are obtained by numerically solving equations (|13p 
and (|14p with values of ci and d\ determined using equa- 
tions ([To]) and (fTE) . Figure [3] indicates that for values 
of (p r ,a r ) same as the initial conditions of figure [3] and 
for (pi = 0.4, a i = 0.01), the density profiles must have 
a boundary layer part near x = 1 merging to constant 
bulk profiles with pi, = 0.4 and 07, = 0.01. This bulk 
part satisfies the boundary condition at x = 0. Numer- 
ical solutions of the hydrodynamic equations for these 
boundary conditions are shown in figures |4] and Fig- 
ure [3] shows that the trajectory with initial condition 



FIG. 6: Phase-plane diagram on the p — a plane for different 
values of m. Trajectories are drawn with initial condition 
p(x = 0) = 0.2 and aix = 0) = 0.8. The arrows indicate the 
direction of more negative value of x. Values of m for various 
lines as we move towards the outermost line are 0.45, 0.58, 
0.64, 0.70, 0.72, 0.73 and 0.74. The flow of the trajectories 
changes at m w 0.725. 



this figure, it appears that there exists a special value 
of m, say, m c , for which the trajectory is a separatrix 
which approaches the saddle fixed point. Our numeri- 
cal solutions show m c w .725. Hence, for m < m c and 
(p r = 0.2, a r = 0.8), the bulk-densities are pb = 0.4 and 
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at, = 0.01. At m = m c , the bulk density value discontinu- 
ously changes to pb ~ 0.597, a h = 0.0125. Naturally, now 
the boundary layers are present at both the boundaries 
to satisfy respective boundary conditions. These bound- 
ary layers are the two separatrices that approach(emerge) 
to (from) the saddle fixed point. Density profiles with dif- 
ferent values of m are presented in figure [7] The inset in 
figure [JJ shows that the boundary layer near x — be- 
comes sharper as the value of e is reduced. The fact that, 
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FIG. 7: Density profiles for different values of m with bound- 
ary conditions (pj — 0.4,(7; = 0.01) and (p r — 0.2, ay = 0.8). 
Values of e for the density profiles with m — 0.72 and 0.73 
are .0.0025 and 0.012, respectively. The rest of the profiles 
are obtained with e = 0.0125. Inset shows the change in the 
shape of the density profile for m = 0.73 as e is changed. Val- 
ues of e as one moves towards the outermost curve are 0.016, 
0.014 and 0.012. 





FIG. 8: Phase-plane diagram for the boundary layer equa- 
tions with non-linear regularization terms. Values of m as 
one approaches the outermost line are 0.63, 0.64, 0.645, 0.648 
and 0.65. Trajectories are obtained with initial conditions 
p(x — 0) = 0.2 and a(x = 0) = 0.8. Arrows on the line 
indicate direction along which x becomes more negative. A 
zoomed version of the same figure is shown in the inset. 



[8j where, from numerical solutions of the hydrodynamic 
equations §§§ and (TTU)) , it has been concluded that in- 
creasing the value of m does not cause a sudden change 
in the value of the bulk density. The present analysis 
brings out the mathematical reason behind such a dis- 
continuous change. Our numerical solutions of the hy- 
drodynamic equations obtained using MATLAB appear 
to be consistent with this observation. 



for m > m c , the bulk density continues to be described 
by a saddle fixed point is clear from the boundary layers 
at both the ends. 

Boundary layers for the equations in © and (fit)]) sat- 
isfy the following first order equations, 

(l-ma)% = p{l-p){l-ma)+c 2 (17) 
ax 

(1 - mp)^- = a(l - er)(l - mp) + d 2 , (18) 

ax 

where, as before, c 2 and d 2 are two integration constants. 
These nonlinear terms do not disturb the boundary layer 
fixed points and since < m < 1 and < p, cr < 1, the 
stability properties of the fixed points also remain the 
same as before. These new factors only introduce minor, 
quantitative changes in the flow trajectories and in the 
value of m c . The flow trajectories for these differential 
equations with the initial condition, (p, cr) = (0.2,0.8), 
and different values of m are shown in figure [8] As a 
consequence, here also the density profiles indicate a dis- 
continuous change in the bulk density as the value of 
m is increased. This result contradicts earlier results of 



IV. MODEL-2: TWO-LANE PROCESS WITH 
PARTICLE EXCHANGE 

Here we consider a system of two lanes with hopping 
of particles from one lane to the other [l2|, EH . Particle 
dynamics on the lanes are governed by the following rules 
(See figure (HJ)). (a) On one of the lanes, say lane-E, 



Lane-Dl 



Lane-E \ 



FIG. 9: Particle dynamics on the two lanes for model-2. 



particles undergo asymmetric simple exclusion process. 
This implies that particles hop to the neighboring site in 
a specific direction provided it is empty. We assume the 
hopping rate to be p here. 



8 



(b) On the other lane, say lane-D, particles are not 
subjected to exclusion interaction and have biased diffu- 
sion with hopping rates to left and right being D~ and 
D + respectively. 

(c) A particle at the ith site on lane-E can hop to the 
ith site of lane-D at rate Ud and the reverse process i.e. 
an attachment of a particle to lane-E from lane-D can 
take place at rate, u a , provided the target site on lane- 
E is empty 

In the steady state, the continuum limit of the statis- 
tically averaged master equations is 



eD 



a dx 2 



8J a 

dx 



Dt - Aa(l - t) = 0, 



(19) 
(20) 



where D T = p/2, D a = (D+ + D~)/2, D = u d N, A = 
u a N and e = 1/N. Here r(x) and a(x) are the average 
densities of particles on the lanes with ASEP and biased 
diffusion respectively with currents on these lanes being 



J T = pr(l — t), and J a 



(21) 



Here, v = D + — D~ denotes the net average velocity 
of particles along lane-D. In addition, we assume that 
the particle reservoirs impose the boundary conditions 
(r(x = 0) = Ti,a(x = 0) = 07) and (t(x = 1) = T r ,a(x = 

l)=<Jr). 

Constant bulk profiles are expected if the terms due to 
particle exchange between the lanes disappear altogether 
jl2l [l3j . This happens under the condition 



D 



A 1 



(22) 



In the following subsection, we consider this model of 
constant bulk profile. Although the phase-diagram for 
the constant profile case has been obtained earlier [ljj, 
this analysis allows us to see how the maximal and min- 
imal current phases and upward and downward shocks 
appear naturally due to the flow of the boundary layer 
solution toward specific fixed points. This study also pro- 
vides an ideal platform to compare this case with that of 
a non-constant bulk profile considered in the next sub- 
section. Therefore, apart from examples with different 
boundary conditions, the next subsection, also contains 
a part in which major differences from the constant bulk 
profile case are discussed. 



A. Constant bulk profile 



Adding equations (fT9"|) and ([2"U]) . and using 
obtain the following steady state equation for r 



, we 



d , dr d 
ox ox ox \ 1 — r 



d_ 

dx 



\vD 



ad 



(i 



+ pr(l-r)]=0, 



(23) 



which should be solved in the presence of the bound- 
ary conditions, t(x = 0) = t; and t(x = 1) = r r . 
Here, D aa \ = D/A. Terms within the second square 
bracket of equation (|23p together give the total current 
•hot = Jt + Ja on t ne t w0 lanes. Special values of t, tm 
and T m , corresponding to maximum and minimum of Jtot 
(see figure fT0|) play an important role in the phase-plane 
analysis. 




"M 



FIG. 10: Jtot is plotted with r. For this plot, 
£>acl = 0.01. 



1 and 



Once again, the boundary layers are the solutions of 
the second order differential equation (|23|) which after 
one integration appears as 



D a D 



pr(l - t)] = 



1 



ad " 



- \vD. 



a< 1 " 



(24) 



Sat- 



Here c is the integration constant and x = — '■ 
uration of the boundary layer to the bulk density, T(,, 
requires 



-vD 



Tb 



ad 



1-71 



pn(l - n). 



(25) 



Since < < 1, the value of c lies in the range — oo < 
c < 0. For a single differential equation as (|2"4"]l . one can 
have a two-dimensional fixed point diagram as shown in 
figure [TT] This is an additional advantage of the present 
system over the previous one. 

The basic principle for finding the density profiles using 
the fixed point diagram is that the shocks or the bound- 
ary layers in the profiles correspond to vertical straight 
lines on the fixed point diagram moving from one fixed 
point branch to the other following the stability proper- 
ties. For example, an upward shock in the bulk of the 
density profile may connect a low-density at the lower 
fixed point branch to a higher density in the middle fixed 
point branch. 
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•0.6 C -0- 4 



FIG. 11: Fixed point diagram on the c — r plane for the 
parameter values v = 1, p = 1 and D a d = 0.01. The arrowed 
vertical lines show the stability properties of various fixed 
point branches. Three fixed point branches, ef, fg and gh 
are referred in the text as lower, middle and upper branch 
respectively. 



(a) Given this rule, the density profile can have local- 
ized shocks under the following conditions, (i) r m o < 
n < t m , t m <r r <r m and c(n) = c(r r ). 

Pictorially, an upward shock is possible if 77 located in 
the lower fixed point branch can be connected to r r in the 
middle branch through only a vertical straight line. This 
constraint is stated mathematically through the above 
equality satisfied by c. Let us assume that for a given 
r r in the middle branch, the localized shock appears if 
T i = T ic (See figure HH). Now for 77 = t;i < n c , the flow 
behavior allows a density profile with a bulk part satis- 
fying the boundary condition at x = and a boundary 
layer near x = 1. This kind of boundary layer is repre- 
sented by vertical line (1) in figure Q21 In this case, it is 
not possible to have a boundary layer at x = since that 
would imply a vertical line passing through 77 and going 
further upward after crossing r r . Although such a line 
passes through both 77 and 7>, it cannot represent a bulk 
profile. Similarly, for 77 = 772 > t\ c , a profile, with bulk 
satisfying the boundary condition at x = 1 and a bound- 
ary layer near x = represented by line (2) in figure [T2l 
is the only solution. 

(ii) A downward shock in the bulk profile can be 
present when r m < 17 < tmo and tm < r r < r rn and 
c(ji) = c(r r ). In this case, a bulk profile of value 77 
at left is joined to a bulk density of value r r through a 
downward shock. 

(iii) Another very special line on which localized up- 
ward shocks can be present is given by i> > r m and 
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FIG. 12: Fixed point diagram of figure QT] Vertical lines are 
numbered in order to refer these in the text in the context of 
different boundary layers and shocks. 



V = T m Q. The shock here joins 77 to the minimum cur- 
rent density r m . The density profile at x — starts with 
a constant value 77 and then reaches another constant 
part T m through a shock, represented by the vertical line 
(3) in figure [T^l and finally satisfies the boundary con- 
dition r r through a boundary layer, represented by the 
vertical line (4) in the figure. 

Approach of a shock or a boundary layer to the bulk 
density is governed by equation (|M|) . The variation of 
a small perturbation St near the saturation to a bulk 
density 77, is given by 

^fin) = [7^^ + P(l - 2r b )]Sr, (26) 

dx (1 - T b y 



where /(r) = D T + D a D a 



Thus the approach 



of a boundary layer to the bulk is in general exponential 
except at special bulk values 77, = r m , tm- The length 
scale associated with the exponential approach to the 
bulk diverges at these special values [U, Gil • 

(b) Boundary conditions for which a minimal current 
phase may appear: The flow lines clearly show that for 
T r > T m and T m o < tj < T m , a shock connecting 77 to the 
middle branch is not helpful for meeting the boundary 
condition at x = 1. In this case the only option is to 
have a boundary layer at x = connecting 77 to r m and 
then a constant bulk profile of density r m followed by a 
boundary layer at x = 1 . Such a boundary layer must be 
represented by line (3) in figure [T2l Note that this is the 
only possibility since any other parallel vertical line with 
c > c(r m ) will correspond to an unphysical solution. This 
explains how a minimal current density profile becomes 
an obvious solution for boundary densities in this region. 
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For the opposite case, t; < r m o and same r r , the density 
should have a constant value 77 with a boundary layer at 
x = 1 represented by line (5) in figure H2l 

(c) Boundary conditions for which a maximal current 
phase may appear: For tm < r; < tmo and r r < tm, 
the only possibility is to have a boundary layer at x — 
represented by line (6) in figure [TSJ This line meets the 
fixed point curve at tm which continues as the bulk den- 
sity. The boundary condition at x = 1 is met by another 
downward vertical line (6'). It is clear that if t; = tmq, 
the density profile can have a localized downward shock 
of height Tj\/ — tm ■ The bulk density will no longer have 
a value tm if U > tmo- In that case, the only possible 
shape for the density profile is a flat profile of density 
T; followed by a downward boundary layer satisfying the 
boundary condition r r . 



B. Non-constant bulk profile 

In this section, we illustrate how a similar analysis as 
that of the previous sub-section can be extended to a 
case of a non-constant bulk profile. In order to produce 
a non-constant bulk profile, we artificially add a term to 
equation (|23p . This kind of term has appeared earlier in 
the steady-state equation ASEP with Langmuir kinetics 
[lf3| . However, for this model, this term only has a math- 
ematical implication of producing a non-constant bulk 
profile. The final steady state equation for the density r 



d , dr „ „ d 
t—[D T — + D a D &d 



dx 
pr(l 



dx 
r)} 



dx V 1 — t 
2t) = 0. 



9 r n T 



Ignoring the second order derivative term in (|23[) , one 
may consider a simplified equation, valid at O(e ), as 



vD 



ad 



+ P (1-2t)]^ + Q{2t-1 
ox 



0. 



(28) 



As mentioned earlier, the solution of this equation de- 
scribes the bulk profile. In principle, one can solve this 
equation explicitly to find the bulk parts of the density 
profile. However, as we shall show, the shape of the pro- 
file can be predicted without explicitly solving this or the 
boundary layer equation. 

After a brief comparison between the two cases with 
constant and non-constant bulk profiles, we discuss a 
few examples with different boundary conditions. To 
show that our approach gives the right profile, we il- 
lustrate these examples with typical density profiles ob- 
tained by solving the hydrodynamic equation numeri- 
cally. We hope that these ideas can be implemented more 
generally for all other boundary conditions that are not 
discussed here. 

Comparison between the cases with constant 
and non-constant bulk profiles: 



The major difference from the previous analysis is that 
for a flat bulk profile, there is no variation in the density, 
once the boundary layer merges to a fixed point. The 
density in the bulk remains constant at this fixed point 
value. This need not be the case here and, in general, the 
density varies along various fixed point branches after 
the boundary layer merges to a fixed point. Thus the 
entire density profile can be predicted, by knowing only 
the slope of the density profile in the bulk and following 
the vertical arrowed lines for the boundary layer parts. 
The slope of the bulk profile can be determined from 
(f2"8"|) and the arrows on the fixed point branches in figure 
1131 represent increasing or decreasing nature of the bulk 
density with x. Since for a flat profile, the density at the 
bulk has to remain constant at 77, t t or at very special 
values like, say, rj\/, the possibility of seeing localized 
shock here is restricted. The shock appears only when 
these values are such that they can be connected through 
vertical lines in the interior of a lane. Such limitations are 
not present here since the varying bulk density may reach 
certain values which can be connected by vertical lines in 
the bulk. It is because of this that seeing localized shocks 
over a region in the phase diagram becomes more likely, 
in general, whenever the bulk density is not constant and 
its slope is appropriate for supporting a shock. Such a 
phase diagram with a wide region of localized shock was 
reported initially in [l6j . 



(27) 




FIG. 13: Fixed point diagram of figure [TT] Vertical lines are 
numbered in order to refer these in the text in the context of 
different shocks and boundary layers. 



n < Tmo and r M < T r < r m : 

There cannot be a boundary layer at x = since that 
would correspond to a vertical flow line passing through 
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ti and the density profile will not correspond to a physi- 
cally meaningful solution. However, a bulk profile satis- 
fying the boundary condition at x — is possible. This is 
an increasing profile along the lower fixed point branch. 
From this lower branch, there are two ways to satisfy 
the boundary condition at x = 1: (a) through a shock 
(represented by a vertical line in the phase plane) con- 
necting the low-density bulk part in the lower branch to a 
high-density bulk part in the middle branch. This latter 
bulk part, with a positive slope (see the upward arrow 
along the middle fixed point branch) finally satisfies the 
boundary condition at x = 1. (b) through a boundary 
layer (represented by a vertical line in the phase plane) 
that connects the bulk density in the lower branch to a 
density in the middle branch. The boundary layer satis- 
fies the boundary condition at x = 1 before it saturates 
to the other fixed point. Figure [14] shows a density profile 
following option (b). 
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FIG. 14: Density profile for n = 0.2 and TV = 0.84. Values 
of other parameters are e — 0.002, v = p = 1, D a d = 0.01, 
= 0.01 and D a = 1. 



Tmo < n < Tm and r M < r r < r m : 

An increasing bulk density satisfying the boundary 
condition at x = may be connected to another bulk 
part in the middle branch through a shock represented 
by the vertical line (1) in figure [TUl This is similar to op- 
tion (a) of the previous example. Such a density profile 
appears as shown in figure [T51 There are still two other 
options each with an upward boundary layer. One pos- 
sibility is that an upward boundary layer at x = meets 
the middle fixed point branch and a bulk part continu- 
ing from there finally satisfies the boundary condition at 
x = 1. The other possibility is that a bulk part with pos- 
itive slope satisfies the boundary condition at x = and 



an upward boundary layer following the bulk satisfies the 
boundary condition at x = 1. 
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FIG. 15: Density profile with a shock for t; = 0.256, r r = 0.81 
and e = 0.015. Values of other parameters are same as those 
of figure [Til 



r r < t m and r M < ri < r m : 

This is a very unique situation for the following reasons. 
Figure IT31 shows that there is only one route on the fixed 
point diagram that the density profile must follow to sat- 
isfy both the boundary conditions. At x — 0, there must 
be a boundary layer of negative slope (represented by 
a vertical line (2)) merging to density tm- This should 
be followed by another downward profile (represented by 
(2')) that satisfies the boundary condition at x = 1. The 
entire density profile is thus described by only two verti- 
cal lines (2) and (2'). A 'bulk' like part of almost constant 
value (~ tm) docs appear in the density profile (see fig- 
ure I16p . However, this is not a real bulk profile since 
tm is not a solution of (|28p . It looks like a bulk density 
since the slopes of the vertical lines (2) and (2') vanish at 
t = tm and hence the profile looks almost flat, r = tm is 
in fact a point of inflection of the density profile. In this 
sense, there is a subtle difference between the maximal 
current phase seen here and in the simple ASEP model 
for which the bulk-density value, 1/2, is a solution of the 
corresponding outer equation [19 ] . 

T\ , T r < r m0 and ti > r r : 

The bulk profile must satisfy the boundary condition at 
x = and increase along the lower branch and then sat- 
isfy the other boundary condition through a boundary 
layer of negative slope at x = 1 , represented by a vertical 
line (3) in the figure [T3l A typical density profile appears 
as shown in figure [T71 
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FIG. 16: Density profile for n = 0.8, r r = 0.1 and e = 0.005. 
Values of other parameters are same as those of figure IT4l 
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FIG. 17: Density profile for n = 0.18, r r = 0.04 and e = 
0.002. Values of other parameters are same as those of figure 

us 



V. CONCLUSION 

We have shown how the phase-plane behavior of the 
boundary-layer equation can be used to understand the 
shape of the steady-state density profile of multi-lane 
driven exclusion processes. We have considered two dis- 
tinct multi-lane processes. In the first model, the lanes do 
not exchange particles but the particle dynamics on the 
two lanes are coupled since hopping of particles on one 



lane is affected by the particle occupancies on the neigh- 
boring lane. In the hydrodynamic approach, the bound- 
ary layers of the density profiles on the two lanes are 
described by coupled nonlinear equations. Using phase 
plane trajectories, we show that, as the interaction be- 
tween the lanes is increased, the bulk density, for certain 
boundary conditions, may increase discontinuously due 
to the influence of a saddle fixed point. This has been 
done using linear phenomenological regularization terms 
for the differential equations. Since this method does not 
require explicit solution of the density profile, it can be 
easily extended to other system of equations with non- 
linear regularizing terms. One such case with nonlinear 
regularizing terms has been discussed here. The second 
model involves two lanes which can mutually exchange 
particles. In addition, particles in one lane go through a 
driven exclusion process and those on the other lane have 
biased diffusion without any exclusion constraint. We 
consider two variants of this model; one with a constant 
bulk density profile and another with a non-constant bulk 
profile. The fixed point diagram shows saddle node bi- 
furcations of the fixed points of the boundary-layer equa- 
tion. For constant bulk profile, the fixed point diagram 
allows us to predict the density profiles and the location 
of the boundary layers uniquely under various boundary 
conditions. For non-constant bulk-profile, one requires 
an additional information regarding the slope of the bulk 
solution. Some of the useful features of the method are 
as follows. With the required information about the bulk 
profile, one can exactly predict under which condition 
boundary layer or shock will have positive or negative 
slope. Similar analysis can be done with ease for different 
parameter values which may reverse the slopes of the bulk 
solutions on some of the fixed point branches or change 
the fixed-point diagram significantly. Analysing the two 
cases of constant and non-constant profiles, it becomes 
transparent as why localized shocks appear only under 
restricted boundary conditions for a flat profile but over 
a wide range of boundary conditions for a non-constant 
profile. 

We have analysed the first model for only some specific 
boundary conditions and we have not made any attempt 
to explore the phase diagram although the method can 
be used in a similar manner for other boundary condi- 
tions also. There are several aspects which may make the 
analysis interesting. One, for example, is the presence of 
more number of physically acceptable fixed points. This, 
apart from leading to interesting phase portraits, may re- 
veal the role of various eigenvalues in deciding the shape 
of the density profile. Second, is the presence of imagi- 
nary eigenvalues which may be responsible for oscillatory 
density profiles. Whether the multi-lane problem consid- 
ered here has this kind of complexity can be explored 
through a detailed analysis for the full phase-diagram. 
This work is in progress. 

For the second model, the condition relating the two 
densities on the two lanes simplifies the problem in two 
ways. First, it leads to constant density profile and sec- 
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ond, due to this relation, it becomes sufficient to con- 
sider hydrodynamic equation for a single density variable. 
In the absence of such a condition, individual boundary 
layer equations look apparently decoupled at the lowest 



order. The boundary layer analysis for the two densities 
would, however, be coupled due to the coupled nature of 
the bulk equations. This may give rise to an interesting 
boundary layer analysis not encountered earlier. 
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